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Single particle dynamics in electron microscopes, ion or electron lithographic instruments, particle accelerators, and particle 
spectrographs is described by weakly nonlinear ordinary differential equations. Therefore, the linear part of the equation of 
motion is usually solved and the nonlinear effects are then found in successive order by iteration methods. 

When synchrotron radiation is not important, the equation can be derived from a Hamiltonian or a Lagrangian. The Hamil- 
tonian nature can lead to simplified computations of particle transport through an optical device when a suitable computational 
method is used. H. Rose and his school have contributed to these techniques by developing and intensively using the eikonal 
method |l[||]. Many ingenious microscopic and lithographic devices were found by Rose and his group due to the simple 
structure of this method 0-^). 

The particle optical eikonal method is either derived by propagating the electron wave or by the principle of Maupertuis for 
time independent fields. Maybe because of the time dependent fields which are often required, in the area of accelerator physics 
the eikonal method has never become popular, although Lagrange methods had been used sometimes already in early days Q. 
In this area classical Hamilitonian dynamics is usually used to compute nonlinear particle motion. Here I will therefore derive 
the eikonal method from a Hamiltonian quite familiar to the accelerator physics community. 

With the event of high energy polarized electron beams || and plans for high energy proton beams || , nonlinear effects in 
spin motion have become important in high energy accelerators. I will introduce a successive approximation for the nonlinear 
effects in the coupled spin and orbit motion of charged particles which resembles some of the simplifications resulting from the 
eikonal method for the pure orbit motion. 

of the 3x1 matrix f. The Jacobian matrix of the function 

^3 



I. INTRODUCTION 

The well known Lagrange variational principle requires 
5 f Cdt = 5 f[p-'q-H]dt = (1.1) 



r(q) is also r since f = J2i=i 4i^qi r = LQ- 

The generalized momentum is p = d~£ — r T (mjf + 
eA) and the variational principle can thus be written as 



with the Lagrangian C, Hamiltonian 7i, and generalized 
momenta p and coordinates q. 

In this principle all variations of q(t) are allowed and 
therefore the Euler-Lagrange equations of motion hold, 

j t dj£(q, 1 1) = %C(«f, 'q, t) . (1.2) 
For relativistic single particle motion the Lagrangian is 



5 J [p r- 1 f-H]dt = 5 J[{m 1 ; r + eA) T rr- 1 'r-'H}dt 

= 5 J [mjv 2 + eA?'r- H]dt . (1.4) 

If only variations 5h=e are considered which keep the 
total energy H. — E constant, the variational principle 
becomes 

K=e [\p-q- K]dt = 5 n=E J p-dq 

[mjv 2 + eP"r\dt = . (1.5) 



C = — nicy c 2 — r + ef ■ A — e$ 



(1.3) 



The variational principle for constant total energy is 
calle d th e principle of Maupertuis. However, in equa- 
tion (1.5) it does not lead to Euler-Lagrange equations 



where the position r(q) is a function of the generalized 
coordinates q. The Jacobian matrix r of this function can 
be written in the form r — (d^r F ) T and has the elements 
d q Ti . In this efficient notation is the transpose 



' ij — ^qj 
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of motion, since not all variations are allowed. 

A particle optical device usually has an optical axis 
or some design curve along which a particle beam should 
travel. This design curve R(l) is parameterized by a vari- 
able I and the position of a particle in the vicinity of the 
design curve has coordinates x and y along the unit vec- 
tors e x and e y in a plane perpendicular this curve. This 
coordinate system is shown in figure [l| The third coordi- 
nate vector e*i = dR/dl is tangential to the design curve 
and the curvature vector is K = —dei/dl. 
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The unit vectors e x and e y in the usual Frenet-Serret 
comoving coordinate system rotate with the torsion of 
the design curve. If this rotation is wound back, the 
equations of motion do not contain the torsion of the 
design curve. The position and the velocity are 

f = xe x + ye y + R(l) , r = xe x + ye y + hlei , (1.6) 

with h = 1 + xk x + yn y . This method is described in 
H and JlO| and is mentioned here since design curves 
with torsion are becoming important when considering 
particle motion in helical wigglers, undulators, and wave- 
length shifters and for polarized particle motion in 
helical dipole Siberian Snakes fH| . 




FIG. 1. Curvatures k x , /% 
ized coordinates x, y, and I. 



of the design curve and general- 



The variational principle ( |l.5| ) for the three general- 
ized coordinates x(t), y(t), and l(t) can now be writ- 
ten for the two generalized coordinates x(l) and y(l). 
This has the following two advantages: a) The parti- 
cle trajectory along the design curve is usually more im- 
portant than the particle position at a time t, and b) 
Whereas 5 n=E does not allow for all variations of the 
three coordinates, the total energy can be conserved for 
all variations of the two coordinates x and y by choos- 
ing for each position r the appropriate momentum with 
rwyv — \f(E — e$(r)) 2 /c 2 — (mc) 2 . We obtain from 
equation ( |Lq ) 



■it 2 d-t 
p ■ dq — o J [■mjv — 



-* (IT 

evA-—}dl = (1.7) 
dl 



Since all variations are allowed, the integrand is a very 
simple new Lagrangian 



L = TwyvJ x' 2 + y' 2 + h 2 + e(x A x + y A y + hA t ) (1. 



which leads to Euler-Lagrange equations of motion 



p x = d x iL , p' x = d x L , 
p y =d y ,L , p' 



d y L 



(1.9) 
(1.10) 



The integral J ' L(j)dl is called the eikonal. 

Since the Hamiltonian formulation is very common in 
the area of accelerator physics, we will show how the 
eikonal can be derived from a Hamiltonian formulation. 

The equations of motion for the three generalized co- 
ordinates x(t), y(t), and l(t) can be obtained from the 
Hamiltonian 



(1.11) 



^Jm 2 c 2 + (p x - eA x ) 2 + (p y - eA y ) 2 + {pi/h - eA{f . 

In the case of time independent fields, TL is the con- 
served total energy E and there are only five indepen- 
dent variables, rather than six. Note that the veloc- 
ity dependent or non holonomic |l3| boundary condition 

W(pt?) Q)i Qi = E cannot be included in the Lagrange 
formalism directly. But in the Hamilton formalism this 
can be done. Furthermore, a switch of independent vari- 
able from t to I can easily be done in the Hamiltonian for- 
mulation. The Lagrange formulation is therefore aban- 
doned (too easily, as will be shown later). In the varia- 
tional condition 



S J [xp x + yp y + Ipi - H]dt = (1.12) 

one can change to the independent variable / as follows: 

S J [x'p x + y% + (-t')H - {-pi)]dl = . (1.13) 

The six canonical coordinates are now x, p x , y, p y , —t, 
and TL, and the new Hamiltonian is given by H = —pi 
which has to be expressed as a function of the six coor- 
dinates 



H 



eAi + 



(1.14) 



(H - e$) 2 - (mc 2 ) 2 - (p x - eA x ) 2 - (p y - eA y f 



In the Hamilton formalism it is simple to take advantage 
of the fact that the total energy is conserved for time 
independent fields; H' = dtH = leads to H = E. Then 
from the six coordinates only the first four have to be 
considered, leading to the Lagrangian 



with df/dl = x'e x + y'e y + hei and dt/dl = \df/dl\/v. 



L = x'p x + y'py - H . 



(1.15) 
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From x 1 = d Po H = ^=(p x - eA x ), y> = d Py H = ^=(p y 
eA y ) where V~ is the square root in H one obtains 



'(£-e$) 2 -(mc 2 ) 2 



= h\ 



X I2 _|_ y l2 _|_ fj2 

h 



77177; 



^^72 + y,2 + h 2 



(1.16) 



Pa; = — x' + eA x , p y = —y' + eAy, (1.17) 



L = m~/vJx /2 + y' 2 + h 2 + e(x'A x + yA y + hA t ) (1.18) 



for 777,77; = {E — e<I>) 2 /c 2 
grangian L agrees with the integrand ( 1 



(mc) 2 . The very simple La- 
of the eikonal. 



In the following it will be shown how the Hamiltonian 
and the Lagrangian equations of motion for the particle 
trajectory q(l) can be solved in an iterative way. We write 
a general equation of motion for a coordinate vector z in 
the form 



(1.19) 



where we assume that z — is a solution of the differen- 
tial equation. Furthermore, we assume z to be small and 
let / be linear in the coordinates. We assume that the 
nonlinear part of the equation of motion can be expanded 
in a Taylor series f^ 2 - The linearized equation of motion 
is solved by a trajectory Zi(l) = M(l)zi which depends 
linearly on the initial coordinates. For the transport ma- 
trix Mil) we therefore have 



M[z, = f 1 Mz, 



(1.20) 



for all coordinate vectors z,; f 1 being the Jacobian matrix 
of/ 



One can write every solution of ( 1.19|) as z(l) 
M(l)^(l), leading to the equation of motion 



M'C + MC' = fMC + f^ 2 (z 



(1.21) 



The Taylor coefficients of Ci^ii with respect to the ini- 
tial coordinates z*i — ((0) are called aberration coeffi- 
cients. With equation (1.20) one obtains 



z(l)=M(l){zi + / M-\l)f^ 2 {z{l))}dl . (1.22) 
Jo 

Now we assume that the general solution z{zi,l) can 
be expanded in a power series with respect to the ini- 
tial coordinates. Then symbolizing the jth order Taylor 
polynomial with [. . we write the orders up to j as 
zj = [z(zi,l)]j, i.e. we use lower indices to describe the 
order of z*i. The upper index in / describes the order in 
z, which is in turn a nonlinear function of Z{. When z n -\ 



is known, one can iterate the expansion up to order n 
with equation ( 1.22| ), since 



z n = M(l){z t + M- x (/)[/^ 2 (4_i)] n }df. (1.23) 
Jo 

The zeroth order of the expansion with respect to the 
coordinates must vanish, which means that the trajec- 
tory q — must satisfy the equation of motion for some 
momentum p(l). Additionally we require that the vector 
potential on the design curve is gauged to zero. This can 
always be achieved. The canonical momentum p then 
also vanishes for the trajectory q — 0. It then follows 
that the Hamiltonian and the Lagrangian have no com- 
ponents linear in the coordinates and momenta. When 
computing trajectories through a particle optical device, 
it is customary to normalize the momenta to the initial 
design momentum pg = p(0). The following two dimen- 
sional generalized coordinates are therefore used: 



x\ fPx/PO 
kVJ ' \Py/Po 

L(q, = L/po , p = d$,L, p' = d$L, 

H{q,p,l)=H/po , q' = d p H : p' 



(1.24) 



The Euler-Lagrange equations lead to the second order 
differential equations ^d^iL — d^L for the two dimen- 
sional vector q. 



II. SUCCESSIVE APPROXIMATION IN TERMS 
OF HAMILTONIANS 

In the Hamilton formalism one obtains first order equa- 
tions of motion for the four dimensional vector z 1 " = 
(<Zi , 32) Pi 1P2) ■ With the antisymmetric matrix J one can 
write the equation of motion as 



-7=1 f 2 ) .:•' = •/'•>•// 



(2.1) 



with the 2x2 identity and zero matrixes 1 2 and 2 . 
This structure implies special symmetries for the trans- 
port maps M of particle optics. These maps describe 
how the final phase space coordinates Zf = M.(z*i) of a 
particle, after flying through an optical device, are re- 
lated to the initial coordinates Z\. These maps are often 
weakly nonlinear and can be expanded in a Taylor expan- 
sion. The Hamiltonian nature implies that the Jacobian 
M. = (dsM. T ) T of any transport map A4(z) is symplectic 
Jl3[ | , meaning that 



M J M = J 



(2.2) 



For the successive approximations we separate the 
equation of motion into its linear and nonlinear part, 



Jdz(H 2 + H^) . 



(2.3) 
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After we have solved for the linear trans port matrix z\ = 
Mzi, we can iterate by equation ( 1.23 ) which takes the 
form 

z n = M{z*i+ [ M-'lJdzH^iz^Udl} (2.4) 



With the relation M 1 J = J M T from equation ( |2~^ ) 
this can be written as 



z n = M{zi + J / [dzH- z {z n -i)] n dl} . (2.5) 
Jo 

The corresponding equation for the aberrations £„ = 
M~ x z n becomes 



Cn 



J / [fcrff^CMCn-l)]** 



(2.6) 



This form of the iteration eq uatio n is quite simple. How- 
ever, since the Hamiltonian (1.14) is a complicated func- 
tion, the evaluation of the four integrals can become very 
cumbersome. 



III. SUCCESSIVE APPROXIMATION IN TERMS 
OF LAGRANGIANS 

In || Rose used a variational principle to derive a suc- 
cessive approximation to nonlinear motion based on the 
eikonal. This method iterates position q and momentum 
p in their nonlinear dependence on the initial position <fj 
and momentum pi . Knowing the order n — 1 dependence 
q n -i and p n -i, one has to compute q^-i by differentia- 
tion of q n -i or by inversion of p = d^:L(q, q' , I). Then 
the eikonal can be evaluated to compute the order n de- 
pendence q n and p n . In general it can be cumbersome to 
compute qJl_i and therefore here we derive a new version 
of the eikonal method, which iterates directly <f n ' rather 
than the momentum. 



In deriving the simple form of equation (2.6), advan- 



tage has only been taken of the symplectic first order 
transfer matrix. We therefore wish to exploit this ad- 
vantage again by working with new coordinates which 
are identical with the canonical q and p up to first order 
so that the new coordinates lead to the same first order 
transport matrix M . To first order one obtains 



P{s) 
Po 



Po 
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+ 2 (q,q') (3.1) 



where p(s) is the momentum of a particle traveling on 
the design curve q = 0, and the upper index 1 specifies 
the part of the vector potential linear in x and y. We 
therefore work with the coordinates 



Po * 



Po \Ah 



(3.2) 



Moreover, it can be shown |l^,||,^7 18 that the contribu- 
tion from the vector potential can be gauged to vanish 
whenever there is no longitudinal magnetic field Bgei on 
the design curve. Then if one investigates trajectories 
which start with momentum po in a region free of such a 
field, we have the simple relation Hi = q,' . 

By splitting the Lagrangian into its second order and 
its higher order part, the equation of motion becomes 



Q' = 



(3.3) 



— — - 1 

p(s) 



p(s) UjJ . _ 



dqL 2 







After having solved the linearized equation of motion, we 
obtain with equation ( 1.23| ) 



Jo \-ji a ?' 
= M{Q l + J J M T (di q 



dl} (3.4) 



T >3 _ Fl-T >3 



An integration by parts leads to 



Q = M{Q % - J 



\M' T ( d «' L 



>3 



M 



-J 



M 







} 



}dl 



(3.5) 



Writing the Jacobian as M T = d^.Qj = d^XQ — Q>2) T 
where Q = Qi + Q>2 was split into parts which depend 
on Qi linearly and nonlinearly, we obtain 

M-'Q = Qi-J [ [(%(?' T - <T>I))%'£- 3 
Jo 



<9 3 (f — qi 2 ))dqL }dl + J 



1 1 



-(d g <g 2 )d ? (L~L 2 )]dl + J 



>:S 



= Qi-l \ [%i- 3 + (%9^)cVL 2 + (%<g 2 K-£ 2 



-J 



{d d £)d r L^ 



(3.6) 



Note that the A' x , A' y , p(s) andpo of e quat ion (3.2) drop 
out of the right hand side of equation (3.6) owing to the 
multiplications by the zeros in equation (3.5). The second 
order L 2 of the Lagrangian is a quadratic form in which 
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every quadratic combination of the q and q' can occur. 
It can be written using a matrix L?_ as L 2 = Q T L 2 Q. 
Part of the above integrand can be rewritten as 

q> T 2 d ? ,L 2 + q% 2 d^L 2 = Ql 2 lZ_Q + Q T LiQ> 2 ■ (3.7) 

For convenience we write L 2 ( a) = a T L 2 a. Another inte- 
gration by parts in equation (3.6) leads to 



J f\d^.(L^ + L 2 (Q> 2 )) 
Jo 



-%Cm 7 £ 2 (Qi) + (%g| 2 )a fl L 2 (Q!)]d/ 



+J 



{d^qU)d r L + {d^q{)d r L^ 



(3.8) 



>3 



The first part of the integral contains Le 
L 2 (Q> 2 ). The integral J Q L E dl is called the perturbation 
eikonal. This scheme embodies the essential requirement 
that the n + 1 order dependence of Le on the initial vari- 
ables Qi can be computed already when Q n -\ is known; 
Q n does not need to be known. For an iteration of Q n , 
knowledge of Q n -i is sufficient. Since q\ satisfies the 
first order equation of motion, we can use the relation 
9q- 1 L 2 (Qi) = ^d^L 2 {Qi) to perform another integra- 
tion by parts, 



M~ X Q 



Qi 



L E dl + J 



(%<?> 2 )%' i 



{d^qi )d ? ,L^ - {d^)d 9 .L\Qx) 



L E dl + J 



(d f)d ? ,L 



-(%gf)%'£ 2 ~ (%gt 2 )%i 2 (0i 



(3.9) 



9q,L e is the part of d^L which up to order n in Q, t 

does not depend on Q n . Similarly the term outside the 
integral is simply the part of {d^Ji r )dqiL which up to 

order n does not depend on Q n . We therefore write 
{d^)d r L - {d^)d r L 2 - (%< 7 | 2 )% i2 ( ( 3i) = 
{{d^(f)d^'L}E- If we now express the Lagrangian in 

terms of the aberrations f with Q ~ M£, we obtain the 
iteration equation 

in =n Qi-Jdn. / L E (M.in-l)dl 



(3.10) 



J 



{{d$f)d r L} E 



When computing £„ from (, n -i with this iteration equa- 
tion, all parts of the right hand side which contribute to 
higher orders are neglected, as indicated by =„. This 
iteration equation can have sev eral a dvantages over the 
Hamiltonian iteration equation (1.23): 



a) The Lagrangian (1.8) is a much simpler function 
than the Hamiltonian ( 1.14 ). 



b) The derivative in equation (3.1C) is performed after 
the integral has been evaluated. Therefore only one 
integral has to be computed and it describes all four 
coordinates of . 

c) The fact that the various coordinates are the 
derivatives with respect to initial conditions yields 
very simple relations |l^] between the various ex- 
pansion coefficients of £„, which are the so called 
aberration coefficients of particle optical devices. 
These relations can be much simpler than relations 
entailed by the symplectic symmetry implicit in the 
Hamiltonian formulation. 



can 



d) The second pair of coordinates in equation (|3 

be calculated very easily. With equations (1.21) 



and (3.3), the equation of motion for £ is 



Af)~, r>3 



%L>3 



(3.11) 



After having computed q n = M 2X 4^ n by iteration, 

the derivative q ' can then easily be computed as 
ti 



q n = M 2xi in using equation ( 3.11 ). One thus only 
needs to iterate the two dimensional vector q n and 
not a four dimensional vector z n as in the Hamil- 
tonian iteration procedure. 



IV. SUCCESSIVE APPROXIMATION FOR SPIN 
ORBIT MOTION 

The time variation of a spin s in the rest frame of 
a particle is described by the so called Thomas-BMT 
equation s = fl but x s polpij where 



BAIT — 



7717 

7 



{(aj+l)B x + (l + a)B l{ 
1 



(4.1) 



(3xE(a+- 
c 1 



'7 



;)} 



with the electric field E, the parts of the magnetic field B 
which arc perpendicular (_L) and parallel (||) to the par- 
ticle's velocity, and the anomalous gyro-magnetic factor 

9-2 
a = 2 • 

Changing to the comoving coordinate system of fig- 
ure |l]), we obtain s = S x e x + S y e y + Si&i and s ' = 
\S X Sin x )e x + (S'y SiKy)ey + (S 1 ; + S x k x + SyKy)ei. 

The equation of motion for the vector S of these spin 
components is then 



S' = Q x S , 



(4.2) 



fi = O.BMT- \/x' 2 +y' 2 + h 2 - K X 

V 
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The equations of motion for the phase space vector z and 
the spin S have the form 



z' = f(z,l) , S' = n(z,l) x S 



(4.3) 



The general solutions transporting the coordinates along 
the optical system, starting at the initial values z\, Si, is 
given by the transport map M. and the rotation matrix 
ReSO(3), 



quaternion A, a spin gets propagated from I to I + dl by 
the quaternion with bo = 1 and b = ^ildl. The resulting 
total rotation is given by A + A'dl and we obtain the 
differential equation 



/ 

fi 2 





-n 2 


-n 3 





-n 3 


n 2 


n 3 





-fix 


Q 2 








(4.9) 



z(l)=M(zi,l) , S(l) = R(zi,l)Si 



(4.4) 



In order to find the general solution, one could compute 
the nine coefficients of the rotation matrix by solving the 
differential equation 



(4.5) 



where the vector product was expressed by the totally 
antisymmetric tensor cuu- However, computing the nine 
components of the rotation matrix seems inefficient, since 
a rotation can be represented by three angles. It has 
turned out [^2j to be most efficient to represent the rota- 
tion of spins by the quaternion A which gives the rotation 
transformation in the SU(2) representation as 



A = a 1 - 



(4.6) 



Here 1 is the 2x2 identity matrix and the elements of the 
vector <7 are the three two dimensional Pauli matrixes. 
When a rotation by an angle (j) is performed around the 
unit vector e, the quaternion representation of the ro- 
tation has ao — cos((/)/2) and a = sin(0/2)e. Therefore 
Og+a 2 = 1 and the identity transformation is represented 
by a = 1. 

If a particle traverses an optical element which rotates 
the spin according to the quaternion A and then passes 
through an element which rotates the spin according to 
the quaternion B, the total rotation of the spin is given 
by 

C = cqI — ic • (7 = (6ol — ib ■ a) (aol — ia ■ a) 

= (boag — b ■ a)l — i(6oa + bao + b x a) ■ a . (4.7) 

The concatenation of quaternions can be written in ma- 
trix form as 

c=r:)=iir) . ,4.m 



B = 



This concatenation of two quaternions can be used to 
find a differential equation for the spin rotation. 

While propagating along the design curve by a distance 
dl, spins are rotated by an angle fldl = \Q\dl around 
the vector fL After having been propagated to / by the 



5) 


= R 


fa \ 


1 






v«; 






-bi 


-b 2 


-b 3 


h 


b 


-bz 


h 


b 2 


b 3 


b 




\h 


-b 2 


h 





Writing the matrix as f2 and the vector as A, the spin 
orbit equation of motion has the form 



z' = f(z,l) , A' = n(z,l)A 



(4.10) 



The starting conditions are z(0) = zq, ao — 1, and a = 
0. The quaternion A depends on the initial phase space 
coordinates Zi and can be expanded in a Taylor series 
with respect to these coordinates. In the following we 
want to devise an iteration method for A n , which is the 
Taylor expansion to order n of A. 

The rotation vector 51 is split into its value on the 
design curve and its phase space dependent part as 
(l(z, l) = f2°(Z)+il- 1 (z, I). The spin motion on the design 
curve is given by A' (l) = Q A (l). Similarly to equation 
(1.21), spin aberrations are defined with respect to the 
leading order motion. Small phase space coordinates will 
create a rotation which differs little from j4o(0 and we 
write the phase space dependent rotation as a concatena- 
tion of Aq and the z dependent rotation (1 + 6,6) which 
reduces to the identity for z = by requiring that the 
aberrations 6 and 6 vanish on the design curve. With 



equation (4.S) we obtain 



(4.11) 



The quate rnion A is now inserted in the differential equa- 
tion (4.10) to obtain 



= (0° + ^)^ 



1 + 5 
6 



.12) 



Taking into account the equation on the design curve and 
the fad 
obtain 



the fact that A^ describes the inverse rotation of Aq, we 



i 



Q(z, I) 



1 



(4.13) 



Writing the Taylor expansion to order n in zi one finally 
obtains the iteration equation 



l ^( 1+ s n 5 T) dL (7 )=0 - (4 - 14) 
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This iteration method was used for the spin transport 
in the program SPRINT |^] and was evaluated using 
MATHEMATICA in @. 

In the case of successive approximation in terms of 
the Hamiltonian, the various aberration coefficients were 
related by the symplectic symmetry. With the Lagrange 
formalism the various aberration coefficients were related 
by their being derivatives of a common integral with 
respect to different initial coordinates. In the case of 
the successive approximation for spin motion, the vari- 
ous aberration coefficients in S and 6 are related by the 
relation (1 + 5f + P = 0. 
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